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lytic upper bounds are obtained for the thermalization rates of the non-local 
observables. 
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1 Introduction 



The study of out of equilibrium processes in quantum field theory is an important and 
interesting problem but the available theoretical toolbox in this area is limited, espe- 
cially when it comes to strongly interacting theories. From the experimental point of 
view, cold atom systems provide a unique setting for exploring the quantum dynamics 
of strongly correlated quantum systems [Tj. As an example, it has become possible to 
set up a quantum phase transition between a superfluid and a Mott insulator, with 
the help of optical lattices [2]. One can even study the quantum dynamics of the 
system as the height of the optical lattice is changed suddenly in time, leading to a 
quantum quench from the superfluid phase to the Mott insulator phase [3j. In this 
manner cold atom systems provide a way of experimentally realizing and simulating 
strongly correlated quantum dynamics near quantum critical points. 

Quenches near quantum critical points are particularly interesting from a theoret- 
ical point of view, because the response of the system is expected to be universal and 
to apply to many different physical systems. Indeed there has been a lot of interesting 
theoretical work studying quench dynamics near quantum critical points |IH7] (for a 
review see [8]). The work in jH[5] provided a nice heuristic picture of the dynamics 
of two point correlation functions and entanglement entropy after a quantum quench 
in 1+1 dimensional conformal field theory (CFT). The heuristic picture, which we 
will refer to as the "horizon" effect, starts with the production of a large number of 
" quasiparticles" from the quench, as the ground state of the initial Hamiltonian is a 
highly excited state of the quenched Hamiltonian. Quasiparticles that originate from 
adjacent points are entangled through the initial ground state. As they travel at the 
speed of light after the quench, they lead to correlations in non-local observables, 
that propagate on a light cone. Inside the light cone, the observables average to their 
thermal values while causality forces correlations outside the light cone to take the 
same form as in the initial vacuum. In [5J it was found that for free lattice systems 
the "horizon" was smoothed out in the sense that, while thermalization started at the 
group velocity of the fastest quasiparticle, complete thermalization was only achieved 
at the group velocity of the slowest quasiparticle. 

Holography provides a novel way to study strongly correlated quantum systems 
and has been used in [9HT5] to study quenches in strongly coupled CFT's in 1+1 
as well as 2+1 and 3+1 dimensions. A similar "horizon" effect has been observed 
within the holographic approach. Alternative holographic models of quenches have 
been constructed and studied in [T6|[TT]. 

The powerful methods of CFT that have been used to study quantum quenches 
apply to "relativistic" theories with an underlying conformal symmetry. Generic 
quantum critical points in condensed matter systems are, however, not conformally 
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invariant but instead exhibit anisotropic scaling of the form, 



x,t) -»■ (Ax, \H) 



(1) 



with a non-trivial dynamical critical exponent z > 1. To study holographic duals of 
these more general quantum critical points it was suggested [18J to consider gravita- 
tional theories in a spacetime that asymptotes to the form3 



where the scaling symmetry (x, t) — > (Ax, X z t) is realized as an isometry of the metric 
when combined with the scaling r — > A _1 r. As is customary, we will refer to a 
metric that asymptotes to fl2]) as an asymptotically Lifshitz metric and the aim of this 
paper is to study quantum quenches for z > 1 using holography in asymptotically 
Lifshitz spacetimes. A priori it is not at all clear whether one should see a similar 
picture arise as in the case of z = 1. In particular, the quasiparticle dispersion relation 
(assuming the concept of a quasiparticle makes sense in these theories), u oc k z , 
suggests that thermalization could start infinitely fast but take an infinite time to 
complete. We will see, however, that this is not what happens in the holographic 
model we consider here. Instead, interactions change the picture qualitatively, and 
we recover a "horizon" effect with a characteristic velocity v oc 

A further motivation for studying quenches in asymptotically Lifshitz spacetimes 
is that some holographic models used to study field theories at finite density, give rise 
to geometries that approach the Lifshitz spacetime ([2]) in the infrared limit of small 
r [2U] . Examples of such IR Lifshitz finite density systems are holographic superfluids 
in their ground state [2TJ and "electron stars" J22j[23]. Thus our results may shed 
some light on the non-equilibrium dynamics of holographic finite density systems. 

A quench can be achieved by suddenly changing the values of coupling constants 
in a given field theory. In holography, the coupling constants of the dual field theory 
are related to asymptotic values of bulk fields as r — > oo. A sudden change in the 
values of coupling constants leads to translationally invariant shifts in the fields near 
r — > oo. Since an asymptotically Lifshitz spacetime acts as a gravitational potential 
well, in a manner analogous to AdS spacetime, the field excitations are pulled towards 
r = and quickly accelerate to a speed very close to that of light. Soon after the 
quench, one is left with a sharp shell of energy density starting from r — > oo and 
falling into the bulk at the speed of light. Eventually the shell of energy will form a 
black hole in the bulk, which corresponds to approaching thermal equilibrium in the 
dual field theory. (There are also alternative ways to achieve a quench holographically, 
see HSHnj.) 

A simple model for a shell falling at the speed of light is given by the Vaidya 
spacetime |24p25] . Indeed the work in [26] showed that for sufficiently small amplitude 

1 See also [19] for early work on gravitational backgrounds with anisotropic scaling. 
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quenches of marginal operators, the spacetime is well approximated by the AdS- 
Vaidya metric. The asymptotically Anti-de Sitter Vaidya geometry was used to study 
holographic quenches in [9T fT2jlT4"lll5p27] . In this paper we generalize the Vaidya metric 
to an asymptotically Lifshitz spacetime. We do this in the context of the family of 
models studied in [28||29]. which have Einstein's gravity, a dilaton and some number of 
U(l) gauge fields in the bulk. We then consider non-local probes (correlation functions 
of gauge invariant operators and holographic entanglement entropy) in the quench 
geometry to uncover the characteristic time and length scales of the thermalization 
process. More precisely, we will analyze lengths of geodesies and minimal surface areas 
and assume that the usual holographic proposal [30H32] extends from asymptotically 
AdS spacetime to Lifshitz scaling so that we can interpret the results as entanglement 
entropy in the dual field theory. Our results may also be relevant for work studying 
entanglement entropy in non-relativistic field theory, see [33H35] . 

The outline of the paper is as follows. In Section 2, we briefly review the black 
hole solutions found in [29] and in Section 3, we construct Vaidya type solutions 
that describe a holographic quench. In Section 4, we study equal time two point 
functions in the quench state, using the geodesic approximation. We also derive an 
analytic upper bound for the velocity at which the thermalization of the two point 
function spreads. In Section 5, we study the entanglement entropy in the quench 
state and again derive an upper bound for the velocity at which the thermalization 
of the entanglement entropy can spread. 



4 



2 Static black holes 



We will consider the gravitational theory specified by the action 
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This theory has a family of Lifshitz solutions. First we will consider the case when 
the number of U(l) gauge fields in the bulk is N = 2. The generalization to a 
larger number of spacetime dimensions and to larger number of U(l) gauge fields 
is straightforward. Throughout we will denote Einstein's equations and Maxwell's 
equations as 



E. 



A" 7 



0, M? = D^e^F^) = 0. 



(4) 



First we will review the static Lifshitz solutions obtained in [29]. The three parameter 
family of solution found in [29] is 



ds 2 = -r 2z b{r)dt l + 



dr 2 



r 2 b{r) 



+ r 2 dy? 
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specified by the parameters (p,,p2,m). For (jSJ) to be a solution to the equations of 
motion, the value of the cosmological constant and the parameters Aj must be fixed 

as 

(2 + z)(l + z) 
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2 ' v^T 

The case p2 = m = corresponds to the Lifshitz vacuum solution. Before constructing 
the infalling shell solutions, it is convenient to write the solution ([5]) in an Eddington- 
Finkelstein-like (EF) coordinate system. 

The ingoing null geodesic is easily found from the metric in §5§ as 



dt + 



b(r) 



-dr = 0. 



(7) 



To define the EF coordinate system, we define a new time coordinate through the 
relation 



-2-1 



dv = dt + 



b(r) 



-dr. 
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Next we exchange the t coordinate with the new v coordinate to obtain the metric 

ds 2 = -r 2z b{r)dv 2 + 2dvdrr z ~ 1 + r 2 dx 2 . (9) 

The gauge fields in the solution ((5]) are of the form A® = A^\r)dt. When going to 
the new v coordinate system these become 

A® = A?\r)(dv- r -^dr) , (10) 
b[r) 

so that there is an A r component induced in the EF coordinate system. For fu- 
ture convenience we note that the A r component can be set to vanish with a gauge 
transformation. Thus, the gauge fields in the new coordinate system have the form 
A = A^dv, where A^ are the fields specified in flSJ). 
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3 The infalling shell solutions 



An infalling shell of massless and pressureless charged matter in an asymptotically 
flat spacetime is described by the Vaidya metric [211[25], which is known analytically. 
The Vaidya metric corresponds to a Reissner-Nordstrom black hole, with the free 
parameters corresponding to the mass and the charge promoted into functions of 
the EF time coordinate v. The Vaidya spacetime is sourced by an explicit energy 
momentum tensor for massless null matter T m ^ and a current density j v oc j r 7^ 0. 
Our task is to find a similar solution as the Vaidya solution, which asymptotes to an 
asymptotically Lifshitz spacetime. Again we will have an explicit source in the vv 
component of Einstein's equations and in the r component of Maxwell's equations, 
while all the other Einstein's equations as well as the matter equations will be solved 
without explicit sources. 

Motivated by the conventional Vaidya solution [25] we choose the following ansatz 
for the metric and the matter fields 



ds 2 = -r 2z b(r, v)dv 2 + 2r z 1 dvdr + r 2 dx 2 , 



b(r,v) = 1 - m(v)r-^ + f(v)— —r~ 2 ( 1+z \ 

Az 

A^ = hl (v)^^^^r 2 ^, (11) 

z 

= h 3 (v)^r 2 ^~\ 

where m(v), f{v), hi(v), /12(f), h^{v) are arbitrary functions of v for the moment. Sub- 
stituting the ansatz (ITT]) into Einstein's equations leads to relations between the above 
functions of v. The vr component of the Einstein's equations leads to 

hl = h }/^~\ f = h 2 2 hf^\ (12) 
The xx and yy components of Einstein's equations are 

E xx = E vy = Y^Tr 2 -*^. (13) 

In order to solve them (for z 7^ 1), we must set 

/i 3 (t>) = const. (14) 

Furthermore we can without loss of generality set h 3 = 1. This also leads to h\ = 1 
and 

/ = 4 (15) 
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We are thus left with two arbitrary functions m(v) and h 2 (v). Since h 2 (v) appears 
only in the combination h 2 (v)p 2 , we see that our solution is simply equivalent to 
promoting m and p 2 into arbitrary functions of v 

ds 2 = -r 2z b(r, v)dv 2 + 2r z ' 1 dvdr + r 2 dx 2 , 



b(r,v) = 1 - m(v)r-^ + ^IJL r -2(*+i) 

Az 

am = ^ES.^^ (16) 



Am 



fir 



z 



We will refer to f liZTj) as the Lifshitz-Vaidya solution. The fields in ( ET|) do not yet 
solve all of the Einstein's and Maxwell's equations, but there are the non-vanishing 
components 

E vv = - ^—EEZ P 2(v)p' 2 (v) + (17) 
M; = r^- Z p' 2 {v). (18) 

These non-vanishing components can be identified as sources from charged infalling 
massless matter with a vanishing pressure. The energy momentum tensor of such 
matter has the form T^ u = pu^u v and a current density = p e u^ with = 5^ v . 
Indeed this is of the form we need to solve the Einstein's and Maxwell's equations 
(HZD and ([H]) with 



-^ r - 2 -» rn/M (19) 

iz rZ 



We have a solution for two U(l) gauge fields, but it is trivial to reduce it to the 
case of a single gauge field studied in j2B] by setting p 2 = 0, which makes the 
field vanish identically. We can also find solutions for theories with more U(l) gauge 
fields in [29J by simply promoting the corresponding free charge densities pi, for i > 1, 
and the energy density m into functions of v. Another direction of generalization is 
to consider different values of the spacetime dimension. This generalization seems 
also to work trivially and we have confirmed this for the bulk spacetime dimensions 
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d = 4,5, 6. The explicit solution in the general case is simply [29] 

ds 2 = -r 2z b(r, v)dv 2 + 2r z ~ 1 dvdr + r 2 dx 2 , 

V p 3 (v) 2 v~^ 



b(r, v) = l- m(v)r-( z+d -V + ,7^f, - r ~ 2 ( d +^ 

v ' ; v ; ^ 2(d - 2)(d + z - 4) 



d r A^ = ^/2(d + z-2)(z-l)fiV^r d+z - 3 , (21) 



d r Af = p,(v)^ 2 —r 3 - d - z , (j = 2, N) 
where now d = 4,5,6. This corresponds to the choice of parameters 



A = -l(d + z-2){d + z-3), Ai = -^2^-|, \ = \Jl Z jZ^' ( 22 ) 

where j = 2, JV. It is straightforward to show that such an ansatz indeed solves 
the equations of motion. 



9 



4 2-point correlation functions 



In this section we study 2-point correlation functions in the quench state. For sim- 
plicity, we restrict our attention in what follows to uncharged black brane solutions 
with = and the metric function 

b(r,v) = 1 -m{v)r- d - z+2 . (23) 

This allows us to focus on the differences that arise between our z > 1 solutions and 
the previously studied z = 1 case without the added complication of non-vanishing 
gauge charge. We expect our calculations to carry over to the charged case in a 
straightforward way. Furthermore we will focus on the case d — 4. Higher dimensional 
cases are discussed in Appendix A, where as an example we show the results for d — 5 
and z — 2. 

For the convenience of numerical computations we define a new radial coordinate 
u—l/r (not to be confused with an Eddington-Finkelstein coordinate), in terms of 
which the metric of interest reads 

ds 2 = -u~ 2z b(u, v)dv 2 - 2u~ 1 ~ z dudv + u~ 2 dx 2 , b(u, v ) = 1 - m{v)u 2+z . (24) 

The calculation of the correlation function is performed in the geodesic approxi- 
mation, which becomes more accurate as the scaling dimensions of the corresponding 
operators in the correlation function are increased. The idea of the geodesic approx- 
imation is to perform a saddle point approximation in the path integral over particle 
paths in the bulk (36], to obtain the bulk Feynman propagator. As was discussed 
in [37] one can obtain the boundary theory correlator by pulling the points in the 
bulk Feynman propagator to the boundary and multiplying by appropriate powers of 
the cutoff e. This leads to the boundary theory two point function in the geodesic 
approximation as given by 

(O(x)O(x')} » r 2 V A /^?? ) (25 ) 

where x M (r) is the geodesic with minimal length and A is the mass of the bulk particle, 
or equivalently the scaling dimension of the dual operator. There is no factor of i 
in the exponent in (j25J) because we are specializing to spacelike geodesies, and it 
is convenient to factor out a minus sign from the square root. By symmetry we 
will choose the geodesies to have y =const, where x = (x, y) are coordinates in the 
transverse plane, and furthermore parametrize it with the coordinate x so that 

u = u(x), v = v(x). (26) 

With this ansatz the particle action becomes 

S = A / dxL = A / dx^u- 2 - U - 2z b(u,v)(v'(x)) 2 - 2 U - z - 1 v'(x)u'{x). (27) 
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Since the Lagrangian L does not depend explicitly on x there is a conserved Hamil- 
tonian ^ ^ ^ 

du' U ~^ dv' V u 2 L ^ ^ 

To find the geodesies we need to set up boundary conditions. In this work we will 
compute only equal time correlation functional The geodesic will start from x = —1/2 
and end at x — 1/2. Furthermore there is a turning point at x = around which the 
minimal length geodesic is symmetric. Regularity imposes the following boundary 
conditions at the turning point 

u '(0) = = v'(0). (29) 

We will furthermore denote the turning point as w(0) = u* and v(0) = t>*. With the 
above boundary conditions we can calculate the value of the Hamiltonian 

H=-u~ 1 . (30) 

As the set of two independent equations of motion we will use one of the Euler- 
Lagrange equations and the conserved Hamiltonian (1281) 

uv" + 2v'v! - m 2 - 1 - ^u z+2 d u (u- 2z b(u, v)){v') 2 = 0, (31) 

1 - 2u l - z u'v' - u 2 - 2z b(u, v) (v') 2 - = 0. (32) 

u 

For the purpose of performing numerics it is convenient to transform to dimensionless 
variables as follows 

u = r Q u, v = TqV, x = r x, (33) 

where r is defined in such a way that the position of the horizon of the forming black 
hole is at u — 1. In this way we can relate tq to the temperature at late time thermal 
equilibrium as 

/ 4vrT \ 1/2 

ro= T7T • (34) 



.2 + z, 

In what follows we will work with dimensionless coordinates and drop the tildes. After 
the change of variables to dimensionless coordinates, the equations for the geodesies 
(13T|) and (|32|) stay invariant in their form except that now 

b(v,u) = l-u 2 +*^$. (35) 
m(oo) 



2 For correlation functions with unequal times in the case of a holographic quench for z = 1 
see [13]. 
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4.1 Lifshitz vacuum 



To begin with, we calculate the equal time correlator in the zero temperature Lifshitz 
spacetime, i.e. we set b = 1. Due to time translational invariance, the geodesic 
has to be independent of time t since "momentum" conservation in the t direction 
forbids the geodesic to turn around in time, which again is necessary to satisfy the 
boundary condition of having both ends of the geodesic at the same time. Thus, the 
geodesic has to be independent of time t and we can use dv = dt — u z ~ 1 du = —u z ~ x du. 
Substituting this into (!3T|) leads to 



-£-(««') = -1, (36) 
dx 



which can be easily integrated into 



u = VD + 2Cx - x 2 . (37) 
Furthermore by imposing the boundary conditions u'(0) = and u(l/2) = we get 



u{x) = ^jj-x 2 . (38) 
Using the conserved Hamiltonian we can write the on shell action simply as 

S — A [ dx-^-. (39) 
J u{xy 

Substituting (13B1 into (1351 leads to 

A ' 2 dx 
P/A - x 



S = Al l (4°) 



This integral is divergent near x = 1/2 and must be regulated. We regulate it by 
introducing a cutoff for x at 1/2 — e. In this way the on shell action becomes 

s=A 7o wi^ =2Atanh (1 -T) K - Alog (i)' < 41) 

The cutoff in x can be related to a cutoff e in the holographic coordinate u through 
u{l/2 — e) = e, which leads to the relation e = y/il. In this way the two point correlator 
becomes 

G 2 (l, t) = (0(-l/2, t)0(l/2, t)} » e- 2A e- s = ^, (42) 

which is seen to be independent of the Lifshitz scaling exponent z. This follows simply 
because the spatial part of the Lifshitz spacetime metric is independent of z. 
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4.2 Thermal equilibrium 

The finite temperature equal time two point function can be obtained by solving (1HT|) 
and (132]) with the black hole metric factor 

b(u) = 1 - u 2+z . (43) 



We can solve (l3Tj) and fl32|) numerically with this choice for 6. The result is shown in 
Fig. [TJ The main point is that the correlator behaves as 




Figure 1: Logarithm of the thermal correlator for different values of z. The different 
curves correspond to z — 1, 2, 6, from bottom to up. The red line is a line with slope 
— 1. The figure shows that the thermal correlation function is fairly independent of 
the value of z. 



G 2 (l,t) oce"^, (44) 

for sufficiently large /. The coefficient £ is the thermal correlation length, which is 
related to the temperature as 



4.3 The quench 

Next we can solve fl3T|) and (1321) for the time dependent Lifshitz-Vaidya background 
with 

b(u,v) = 1 -m(v)u 2+z , (46) 

where we choose the profile 

m (v) = ^(1 - tanh(v/v )). (47) 
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For the time scale Vq appearing in (jlTl) we will choose the value Vq = 10~ 2 , or in 
dimensionful coordinates Vq = 10~ 2 (2 + z)/(4irT). In the dual field theory this corre- 
sponds to adding a pulse of energy to the vacuum at the time t = 0. The time scale in 
which the energy pulse appears is vq. We call this operation a quench. Before solving 
the equations numerically we can understand some of the main features without the 
detailed calculation. At times earlier than t = 0, the spacetime will look like the pure 
Lifshitz vacuum and thus, the equal time correlator will take the vacuum form (|42p . 
Also for times much larger than t = the spacetime looks like the Lifshitz black 
hole so that one might think that the equal time correlator will take the thermal 
equilibrium form (jHJ). This conclusion does not generally hold for late times. In- 
stead the qualitative behavior of the correlator depends on the transverse separation 
I. If I is sufficiently small, the geodesic will not "drop" too close to the horizon and 
v(x) > along the entire geodesic. Then the geodesic will indeed be that of the 
late-time black hole spacetime and one recovers the thermal correlator (jHJ). On the 
other hand, when I is sufficiently large, the geodesic passes through both the event 
horizon and the apparent horizon [9], as discussed below, and part of the geodesic 
will have v(x) < 0. This means that the geodesic passes through the infalling shell 
at v = and extends into the part of the spacetime with vacuum geometry. In this 
case the correlator will not be thermal. 




Figure 2: In both of the figures the dot-dashed curve (red) corresponds to the apparent 
horizon, the dashed curve (blue) to the event horizon, and the solid curve (black) to 
a generic geodesic contributing to the two point function. The left hand side figure 
corresponds to z = 2 and the right hand side figure to z — 3. The geodesic is seen to 
pass through both horizons, once through the event horizon and twice through the 
apparent horizon. 

We can conclude that thermalization does not happen globally for the system at 
any finite time. Rather the thermalized region, as seen by the two point function, 
expands in time. This is the "horizon effect" observed in two point correlators in [5j[7]. 
The same effect has been observed before in holography in [T0HT3] for the relativistic 
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case z = 1. The "horizon effect" seems more surprising for other values of z. 

In Appendix B we find the location of both the event horizon and the apparent 
horizon in the infalling shell spacetime with z > 1. In an evolving geometry the 
apparent and event horizons can be very different and this is also seen here. As 
emphasized in [9] , the two point correlation function probes regions of the spacetime 
that are inside both the event horizon and the apparent horizon. Indeed this seems 
necessary to reproduce the "horizon effect" known in quantum field theory. Figured] 
shows examples of geodesies, obtained numerically for different values of z > 1, that 
pass through both horizons and also through the infalling shell into the v < region. 




Figure 3: Lower bounds for thermalization times for different values of z. The data 
points are obtained by numerical integration of the integral in (I49p . From top to 
bottom the different colored data points correspond to the values z = 1,2,3,4,5. 
The lines with the corresponding colors are best fit lines with slopes given by f l54l . 



Thermalization of the two point correlation function spreads with a finite veloc- 
ity. Another way of saying this is that the thermalization time scale t t h erm (l) is an 
increasing function of I with a finite slope. Here we define the thermalization time as 
the time when the two point function reaches the thermal value. First we can obtain 
a lower bound for the thermalization time ttherm(l) in the case of a sharp quench. The 
correlator can be thermal only when the turning point of the geodesic u* is below the 
instantaneous position of the matter shell «o(tjj- We can obtain uo(t) by integrating 
(J7J. This leads to 

'=/. ^wy (48) 

where we use b(u) = 1 — u 2+z , since we are interested in the region outside the matter 
shell. The geodesic equations in the black brane backround can be integrated to 

P* du , s 

1 = 2 (49) 

J° ^(g-l)6(«) 



3 So if the geodesic passes through the shell, the two point function will depend on time, and will 
thus not be thermal. 
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A lower bound for the thermalization time is when uq = since for smaller times 
t (or larger distances /) the geodesies necessarily pass through the shell and make 
the correlator time dependent. Note that this does not necessarily determine the real 
thermalization time as it can be that even at later times there is a geodesic which 
passes through the shell and has a shorter length than the one that probes the region 
outside the shell. The integral (|48p can be seen to be an integral representation of 
the incomplete Beta function 

while we have not found a closed form for the integral in ( )49l) . Still we can find 
analytic expressions for both of the integrals in the limit as t — > oo and I — > oo. 
This limit corresponds to integrating (jlSjl close to the pole at it — > 1, which leads 
to a logarithmic divergence. Similarly the integral in ( H9l) diverges logarithmically as 
it* — > 1. This leads to 

1 

log(l -m ) + finite (51) 



2 + z 

1 ' log(l -«♦)+ finite. (52) 



2 y/2(2 + z. 

Identifying uq = leads to a lower bound for the thermalization time 

t = — . — h finite, (53) 

for large t and I. This gives us an upper bound for the velocity the thermalization 
can spread with 

v = 2^1 + ^/2, (54) 

for large t and I. Transforming back to dimensionfull coordinates using ( |33l) gives 
the dimensionfull velocity (recalling that in a theory with Lifshitz scaling symmetry 
velocity is indeed dimensionfull when 2^1) 



■" = 2 (jTl) V 1 + 2' (55) 

This is one of the main results of this section. Because of the finiteness of the bulk 
velocity of light (or causality in the bulk), there is an upper bound on how fast the 
two point function can thermalize, for all values of z. This upper bound depends 
on the state of the system after quench explicitly through the late time equilibrium 
temperature T according to fl55|) . Also, this velocity is independent of the scaling 
dimension of the operator in the correlation function^ A similar bound for general 
bulk spacetime dimension is calculated in Appendix A. 



4 As long as the scaling dimension is sufficiently large for the geodesic approximation to hold. 
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One should note that the fact that I is linear function of t in (1531) depends only on 
the presence of a first order zero in the metric function b at the horizon. Thus, the 
"horizon" effect (that the thermalization of the two point function spreads in a cone 
with finite velocity) in a quantum quench is indeed directly related to the formation 
of a black hole horizon in the gravitational dual. 

Even though the relation (|53|) was derived in the limit of large I and t it can 
be seen to apply well for sufficiently small I and t as can be seen by evaluating the 
integrals (HHj) and (149]) numerically. This comparison to the numerical evaluation is 
shown in Fig. [3] 




Figure 4: Logarithm of the quench correlator for z=l, with the vacuum value sub- 
tracted. The red line corresponds to twice the speed of light. The blue surface 
corresponds to geodesies passing through the matter shell, while the green surface 
corresponds to geodesies probing the u < 1 region. 




Figure 5: Logarithm of the quench correlator for z=2, with the vacuum value sub- 
tracted. The red line is a reference line with slope dl/dt = 2. 
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Figure 6: a) Logarithm of the quench correlator for z=3 with the vacuum value 
subtracted. The red line is a reference line with slope dl/dt = 2. The thermalization 
is seen to happen very suddenly as there are several geodesies contributing, b) The 
blue solid line is the real thermalization time extracted from the correlator, while the 
red dashed line is a reference line with slope 1 and the brown dot-dashed line is the 
lower bound for the thermalization time as obtained from numerical integration of 
gHD and (P]l. 

Numerical results for the correlation functions for different values of z are shown 
in Fig. HJ Fig. [5] and Fig. [6j For illustrational purposes we plot the logarithm of the 
correlator with the vacuum value substracted. 

To obtain the real thermalization time one has to take into account all the possi- 
ble geodesies and pick the ones that have the lowest length to obtain the correlation 
function. For the ranges of t and I we have studied, the lower bound for the thermal- 
ization time agrees wtth the teal thermalization time for z = lg A s can be seen from 
Fig. [5] and Fig. |6]0 the real thermalization time for z = 2, 3 is indeed bigger than 
the lower bound as there are shorter geodesies available that pass through the shell. 
This is seen in the figures as the lengths of two branches of geodesies crossing. In the 
figure, the surface which is above, has a smaller geodesic length and corresponds to 
the real value of the correlation function. 

It seems that there are two competing effects at work as z is increased. On one 
hand, the upper bound for the thermalization velocity increases as v oc y/l + z/2, 
which can lead to faster thermalization as z is increased. On the other hand it seems 
that for larger z, the real thermalization velocity due to competing geodesies is getting 
smaller than the upper bound leading to a slowing down of the velocity as can be 
seen from the right hand side of Fig. EJ 

5 There is nothing that guarantees this for larger and larger values of I. 

6 It should be noted that Fig. has been made by approximating the time dependent part of 
the correlator (the blue surface in Fig. [6^ ) as being independent of I. This is seen to be a good 
approximation at least for I > 4. A similar approximation has been used to produce Fig. [7] and in 
the next section Fig. |5b and Fig. [TUJ 
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A-hog(G 2 /Gl) A- 1 log(G 2 /G[) 




Figure 7: Logarithm of the two point correlator, with its thermal value subtracted. 
The different curves correspond to different values of I as I = 4, 8, 12 from bottom to 
top. The left figure is for z = 2 while the right figure is for z = 3. 



From Fig. [5] and Fig. [6^, we see that to a good approximation the ratio of the time 
dependent correlator and the vacuum correlator is independent of I (for sufficiently 
large I) at a fixed time to- This tells us that outside the thermalized region the 
correlator can be approximated by 

G a (l,*>)«^, (56) 

where the wavefunction renormalization factor Z is a decreasing function of the time. 
For the cases z = 2,3, the correlation functions at fixed I approach their thermal 
values linearly to a very good approximation as can be seen from Fig. [71 Thus 
we see that the time dependence of the wavefunction renormalization is to a good 
approximation given by 

Z(t„) « e-*» A / T , (57) 

where r is simply the inverse slope of the linear parts of the curves in Fig. [71 Also, 
Fig. [7] shows that the time derivative of the correlator, as it reaches the thermal 
value, is discontinuous. 
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5 Entanglement entropy 



Another probe of correlations in the quench backround is provided by the entangle- 
ment entropy. A holographic formula for the entanglement entropy for a region S 
in the boundary field theory has been suggested in [30H32] to be given by the area 
of a minimal surface that ends on a curve C = dS at the AdS boundary. We will 
use the same prescription to calculate the entanglement entropy in the case of an 
asymptotically Lifshitz spacetime with z ^ 1. On the rest of this section we will 
simply assume that the above prescription for calculating the entanglement entropy 
extends to theories with Lifshitz scaling. 

We will study the entanglement entropy of an infinite strip with a width I. The 
entanglement entropy is now obtained by minimizing the area functional 



Sent = ^J d °\J\ tetd a X*d p Xp\. (58) 

For the infinite strip geometry we will choose to parametrize the surface with coordi- 
nates a 1 = x, a 2 = y where the strip will be infinite in the y direction and has a length 
I in the x direction. Due to translational symmetry of the strip in the y direction, we 
can take an ansatz with d y X^ = 0. The coordinates of the surface are taken as 

u = u(x), v = v(x). (59) 

Substituting this into the area functional gives 

Sent = -^7- I dy I dxu' 1 V ' -{v') 2 u~ 2z b(u, v) - 2v'u'u- z - 1 + u~ 2 . (60) 

4CrjV J J 

It is convenient to define an entanglement entropy density as 

dxu^ 1 \J —{v') 2 u^ 2z b{u,v) — 2v'u'vr z ~ 1 + i ur 2 
dxLA- (61) 



AG]s[S en t 

^ent ~ 



fdy 



The fact that we multiply the entanglement entropy by Gn can be interpreted in the 
dual field theory as dividing by the number of degrees of freedom in the dual field 
theory0 Again there is a conserved Hamiltonian following from the fact that the 
integrand in the area functional does not have explicit dependence on x 

dL A , dL A , 1 

We will choose the x-coordinate in a way that the center of the strip is at x = 0. 
Again the equations and boundary conditions are symmetric under x — > —x. Thus, 



7 In the CFT case, by the central charge. 
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the minimal length geodesic is again symmetric around x = 0. Regularity requires us 
to set tt'(O) = v'(0) = 0. Denoting the turning point of the surface as (u(0),v(0)) = 
(u*,v*), we obtain the value of the hamiltonian (1621) as 

H A = -u-\ (63) 

This together with the Euler-Lagrange equation that follows from extremizing the 
area functional in (JBTT) gives rise to the two independent equations 



uv" + Av'v! - 2u z ' x + (z + l)u l - z (v'f - ^m{v)u 3 {v'f = 0, (64) 
1 - 2m 1 " W - u 2 - 2z b{u, v) {v' f - ^ = 0. (65) 

5.1 Lifshitz vacuum 

In the Lifshitz vacuum 6=1 and we can set dv = —u z ~ 1 du. In this way fl65|) becomes 

1 + (u 1 ) 2 = % (66) 



which can be integrated as 



I [ u * u 2 du ^(3/4) 



2 J v^ 3 ^ * r (V4) 
The entanglement entropy density is now given by 

r 112 ^ u 2 r u * duv~ 2 

s ent = 2 dx^ = 2 (6g) 

where e is defined through u(l/2 — e) = e and we used du/dx = —\Ju*/u A — 1. 
Performing the integral in (1681) we get 

_2 1 ^-1/4)1X3/4) 
^ t - 7 + y f ^ • (69) 



This result is independent of z and has the same value as in AdS spacetime |31j . 
since again the spatial part of the Lifshitz vacuum metric is identical to that of AdS. 
The first term in ( |69l) is divergent as the boundary theory cutoff e is sent to zero. 
Furthermore, since this term comes from the small u behavior of the surface, it is 
independent of the state of the system. Meaning that it has the same value also in 
the thermal state and in the quench state. Thus, from now on we will substract the 
cutoff dependent term out of the entanglement entropy and define 

2 

As ent = s ent - -. (70) 
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The main physical physical point here is that 

As ent OC -y, (71) 

which means that the entanglement entropy of a strip of width I increases as I is 
increased. Since As en t is the entanglement entropy per unit length in the y-direction, 
f JTTj) is indeed consistent with dimensional analysis. 

5.2 Thermal equilibrium 

To study the entanglement entropy at finite temperature we set 

b=l-u 2+z , (72) 

in equations fl64|) and f )65|) . Next one can numerically integrate the equations. Results 
from a numerical integration are shown in Fig. [HJ 




Figure 8: The finite part of the entanglement entropy density as a function of the 
width of the strip I, in the thermal state. The different curves correspond to the 
values z = 2, 4, 6 from top to bottom. The straight line is a reference line with slope 
1. 

The main physical points here are that the entanglement entropy for small dis- 
tances behaves as in the vacuum 

As ent oc - j, I < T- x '\ (73) 

while for large distances the entanglement entropy seems to be a linear function of / 

As ent oc IT 2 I\ I > T- X l\ (74) 

so that it seems extensive at large distances just as the usual thermal entropy. 
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5.3 The quench 

To calculate the entanglement entropy in the quench state we should solve equations 
fSD and (ESD for 

b=l- m(v)u 2+z , (75) 

where m(v) will again be chosen to be a hyperbolic tangent as in ( 1471) . With the 
same arguments as we used in the case of the geodesic calculation in the previous 
section, also the entanglement entropy will thermalize only locally at any fixed value 
of time. Meaning that at large enough values of I the entanglement entropy will be 
non-thermal and time dependent, no matter how large values of t we look at. This 
again leads to a horizon effect similar to that in the two point function. Such a horizon 
effect has been observed before in the entanglement entropy in 1+1 dimensional CFTs 
in [1] and later by using holography in [9HT2]. 

In fact the thermalization of the entanglement entropy will happen slower than 
that of the two point function in the holographic setup. We can see this simply by 
looking at the area functional in ( I6ip . Essentially this differs from the geodesic length 
functional simply by having a one extra power of u~ l in the integrand. This means 
that in order to minimize the area functional, the turning point will want to be at 
a larger value of u than that in the case of a geodesic. Thus, we see that the minimal 
surface will "hang" deeper in the bulk, and will pass through the matter shell with 
a smaller value of I. This means that the thermalization time for the entanglement 
entropy ought to be larger than that in the case of the two point function. 

We can again make this intuition more precise in the case of a sharp quench 
(meaning m(v ) oc 0(v)). The calculation is identical to that performed in the previous 
section except that now we the minimal surface equation integrates to 

1 dU (76) 



2 



3 ~ !)*(«) 



This integral again diverges logarithmically as « t 1 and we can extract the diver- 
gent part as 

2 = ~2^/TTz l ° g(yl ~ U *" > +finite ' ( 7? ) 

Identifying u* with the position of the matter shell given by f l5Tj) . we obtain an upper 
bound for the thermalization time (as the time when the minimal area surface starts 
to penetrate the matter shell) as 

1 I „ 

t = . - + finite. (78) 
l , z 2 

Again this gives us an upper bound for the thermalization velocity (restoring the 
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Figure 9: a) The entanglement entropy density in the quench state as a function of 
the width of the strip I and the time t, for the case of z = 2. The red line corresponds 
to the upper bound velocity, which in dimensionless coordinates is simply 1. b) The 
blue solid curve is the real thermalization time, while the brown dot-dashed curve is 
the lower bound for the thermalization time and the red dashed curve is a reference 
line with slope 1. As the figure shows, thermalization of the entanglement entropy 
really happens slower than the upper bound due to multiple surfaces contributing to 
the same values of (l,t) 



dimensionfull quantities from ( )33l) ) 



The generalization of ( I79p to general bulk dimensionality is shown in Appendix A. 
This velocity is indeed smaller than that obtained for the thermalization of the two 
point function in ( |55|) . Again we should emphasize that (!79|) is only an upper bound 
for the velocity the thermalization can propagate with. To obtain the real thermal- 
ization time one has to do the full numerical calculation of solving equations f !64p and 
fl65|) . The results from a numerical integration are shown in Fig. [9] for the case of 
z — 2. In this case we again see that the real thermalization velocity is smaller than 
the upper bound. 

From Fig. [TD] we see that the entanglement entropy increases to a very good 
approximation linearly with time. Also we see that the entanglement entropy reaches 
its thermal value very abruptly. Meaning that dtS en t is discontinuous right when the 
entanglement entropy reaches a thermal value. A similar effect is seen for the case of 

z = i in dMg. 




(79) 
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Figure 10: Evolution of the entanglement entropy density in the quench state for 
z = 2, with the thermal value subtracted. The different curves correspond to the 
values / = 4, 8, 12. 



6 Discussion 

In this paper we have studied quenches in theories with a scale invariance that is 
anisotropic between space and time. The study of quenches in quantum field theory 
has mainly concentrated on 1+1 dimensional systems near a quantum critical point 
with relativistic scale invariance z = 1, where one can use conformal field theory 
methods [UEJ[7] to study the system. There one finds that after a quench the non- 
local observables such as two point correlation functions and entanglement entropy 
thermalize only locally at any fixed value of time. The size of the thermalized region 
is found to expand with twice the speed of light. This was called the "horizon" effect. 
More generally in "free" lattice systems it was found that the thermalized region 
expands with twice the group velocity of the slowest quasiparticle. This provides a 
nice intuitive picture of how the quench is followed by a release of a large number 
of energetic quasiparticles. As being released from sufficiently close to each other 
the quasiparticles are entangled. Subsequently these quasiparticles move with their 
corresponding group velocities (which in the case of a CFT is the speed of light). 
When a pair of quasiparticles reaches the points where the operators of the two point 
correlation function are located, one finds quantum correlations that are due to their 
initial entanglement. When the pair of quasiparticles has already passed through 
the two points, the correlation function is affected only by non-correlated quasipar- 
ticles, and it averages to a thermal value. When studying quantum quenches using 
holography a similar "horizon" effect was found in 1+1 as well as higher dimensional 
systems [5HT3]. albeit with some small differences. 

The above results are only applicable to the case when the quantum critical point 
exhibits "relativistic" scaling z = 1. In this paper we have studied quenches on 
quantum critical points with general values of z using holography. We constructed 
analytic solutions corresponding to an infalling shell of massless and pressureless 



25 



matter. These solutions are straightforward generalizations of the Vaidya solution [25] 
to an asymptotically Lifshitz spacetime. This solution provides a holographic model 
of a quench where energy is injected to the system at time t = 0. Subsequently the 
system evolves into a thermal state, which in the gravitational description means that 
the matter shell forms a black hole. 

In a quantum critical point with z > 1, one expects the "quasiparticles" (if such a 
concept makes sense) to have a dispersion relation of the form u oc k z . This leads to 
a group velocity v oc k z ~ x . If one applies the heuristic picture of the "horizon" effect 
to these "quasiparticles", one concludes that thermalization of non-local observables 
never happens as the group velocity of the slowest modes is arbitrarily close to zero. 
On the other hand, if the thermalization is somehow driven by the high energy modes 
one might think that thermalization happens immediately as v is not bounded from 
above for large k. 

Neither of these possibilities is what we find in the holographic model. What we 
find is that there is a "horizon" effect which proceeds at a finite velocity v oc 
where T is the temperature of the equilibrium state reached. The peculiar tempera- 
ture dependence of the propagation velocity of the thermalization "horizon" follows 
simply from dimensional analysis. It also suggests that the thermalization is mainly 
mediated by "quasiparticles" with average energy T as this can lead to a velocity 
v oc T^ -1 )/ 2 . This is of course natural since after the quench, the average energy of 
a "quasiparticle" should indeed be given by T. We do not find any time dependence 
in the correlation function or the entanglement entropy after the thermalization time 
is passed. So it seems that there are no "small momentum quasiparticles" that could 
destroy the thermalization, as they would lead to late time oscillations in the correla- 
tors [7]. A possible explanation is that any "quasiparticles" with small momenta will 
be excited to have momenta of the order k = T x l z by scatterings with the surrounding 
medium, as we are dealing with theories which are strongly interacting. Then because 
the "quasiparticle" group velocity is a non-trivial function of the momentum (this is 
different from the case z = 1 where the velocity is a constant, the speed of light) 
the "quasiparticles" can speed up to v oc as given by the average energy 

in the system. Similarly, high energy "quasiparticles" can slow down to velocities 
v oc T^ z ~^' z through scatterings. 

We also found upper bounds for the velocity of the thermalization "horizon" 
analytically. These indeed show that thermalization must happen at a finite velocity 
even though there is no reason for this from causality in the dual field theory for the 
case of z > 1. The velocity bounds, as well as the real thermalization velocities, are 
different for different observables. We found that generally entanglement entropy is 
the observable that thermalizes last. The ratio between the thermalization velocities 
of two point functions and entanglement entropy was found to be order 1. 

Even though the different observables thermalized with different velocities, there 
was some universality in the two point functions because the velocity at which they 
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thermalized was seen to be independent of the scaling dimensions of the operators 
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This may be an artifact of the geodesic approximation. 
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A A higher dimensional example 



In the bulk of the text we were dealing with the case where the bulk was 3 + 1 
dimensional. The results for higher dimensional cases are qualitatively similar. In 
this appendix we consider d dimensional bulk, and present numerical results for the 
case d — 5 and z — 2. 

The geodesic equations (}3"Tj) and (1521) depend on the bulk dimensionality only 
through the metric factor 

b(u,v) = 1 - m(v)u d - 2+z . (80) 

Thus, they generalize immediately to general values of d. By considering the case of 
a thin shell, we can again calculate analytically a lower bound for the thermalization 
time at large values of I. The integrals ( 1481) and (T4TJ1) apply to the higher dimensional 
case with the difference that the metric factor is b = 1 — u d ~ 2+z . With this difference, 
the logarithmically divergent parts are now 

£ = -- ; — - — -log(l-u*)+ finite, (81) 
d + z — 2 

\ = ~ /OM f l0g(1 - M * } + (82) 

2 ^2{d + z - 2) 

transforming back to dimensionfull quantities we get an upper bound velocity (at 
large /) with which the thermalization can spread as 



As an example for the full correlation function the case d = 5 and z = 2 is displayed in 
Fig. [TTJ Thermalization is in this case seen to spread with the upper bound velocity, 
which in this case is (in dimensionless coordinates) v = \/T0. 
The entanglement entropy is now obtained by minimizing 



Sent = dd ~ 2 °\j\ de * 0«A>fyjg. (84) 

Here we calculate the entanglement entropy for a "belt" region. In this way we will 
choose the coordinates parametrizing the hypersurface as a 1 = x, a 2 = y, a d ~ 2 = w. 
The interior of the "belt" corresponds to x G (—1/2,1/2) and the other coordinates 
range from — oo to +oo. In this way the entanglement entropy density functional to 
minimize becomes 

AC S 1 C 

Sent = r , , = / dxu 3 ~ d J-{v') 2 u- 2z b{u, v) - 2v'u' U -^ + u~ 2 . (85) 
J dy...dw J 
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Figure 11: a) The logarithm of the two point correlation function with vacuum value 
subtracted, for z = 2 and 3+1 dimensional dual field theory. The red curve corre- 
sponds to the upper bound speed v = y/lO from fl53"j) . b) The entanglement entropy 
density for the same d and z. The red curve corresponds to the upper bound speed 
v = y/10/3 from (159"]) . 



Minimizing ( )85l) leads to the equations 

uv" + 2(d - 2)v'u' -{d- 2)u z ~ l + {z + d- ^-"{v'f 
+ ^(4-z-d)u d - 1 m(v)(v') 2 = 0, (86) 

1 - 2u l - z v!v' - u 2 - 2z b(u,v)(v') 2 - ^— = 0. (87) 

u 4 

Again we can analytically calculate a lower bound for the thermalization time. This 
time the minimal surface equations lead to an integral 

I r* du 1 , /, x * • 

: log 1 — «*) + finite, 



which together with ( 1ST]) leads to an upper bound velocity (in dimensionfull coordi- 
nates) 

/4vrT\ £ f i ld+z-2 . . 

v - 2 (jti) iw^r (89) 

Curiously, the ratio of the upper bound velocities for the two point function and the 
entanglement entropy is independent of z and is given by yd — 2. The entanglement 
entropy for z = 2 and d = 5 as a function of / and t is shown in Fig. ITTb . In this case 
thermalization happens slower than the upper bound velocity. 



B Apparent horizon and event horizon 

In this appendix we locate the apparent horizon and the event horizon for the infalling 
shell metric 

ds 2 = -r 2z b(r, v)dv 2 + 2r z ~ 1 dvdr + r 2 dx 2 , (90) 
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where we take the number of spacetime dimensions to be d = 4. The calculation 
follows closely the corresponding calculation for z = 1 in [9]. 

Our spacetime can be foliated by the spacelike surfaces r =const and v =const 
and then we can define two null vectors that are orthogonal to the surfaces 

N m = -r l ~ z d r) N out = d v + K z+1 b(r, v)d r , (91) 

where the normalization is chosen in a way that iVj n ■ N out = — 1. The expansion 
around the null directions are given by 9 = g^V^Nu, where g^ u is the induced metric 
on the surface. A straightforward calculation gives 

6 m = -2r- z , 6 out = r z b(r,v). (92) 

Thus, the inward null geodesies always converge, while the outward geodesies converge 
when b(r,v) < 0. The apparent horizon is in this case given by the locus where 
b(r,v) = 0. For the uncharged black hole solutions, this gives the position of the 
apparent horizon as 

r(v) = m(v)2+^. (93) 

Now we calculate the position of the event horizon. Since in our case the mass 
function m(v ) approaches a constant masi)4 oo, the event horizon is simply given 
by the outgoing null geodesic which approaches the position of the apparent horizon 
as v — > oo. Outgoing null geodesies satisfy the first order differential equation 

nv I 

j- = -r z+1 b(r,v). (94) 
dv 2 

For a given mass function the solution to this equation, with the boundary condition 
lim^oo r(v ) = m l ^ 2+z \ specifies the event horizon uniquely. The event horizon and 
the apparent horizon are shown in Fig. [21 for the choice of mass function in f T4T|) . 
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